On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments

10/19/2019

Authors:

Colin Grudzien^1,3, Marc Bocquet² and Alberto Carrassi^3,4

University of Nevada, Reno, Department of Mathematics and Statistics
CEREA, A joint laboratory École des Ponts Paris Tech and EDF R&D, Université Paris-Est, Champs-sur-Marne, France.
Nansen Environmental and Remote Sensing Center, Bergen, Norway
Mathematical Institute, University of Utrecht, the Netherlands

Outline

Intro to data assimilation (DA)
Twin experiments in toy models
The Lorenz-96 model
Multiscale model reduction and random dynamical systems
Stochastic integration schemes and bias in model statistics
Numerical benchmarks:

Strong and weak convergence
Ensemble forecast statistics
Filtering statistics

Conclusions and future work

Data assimilation in the geosciences

Data assimilation (DA) and ensemble-based forecasting are the prevailing modes of:

prediction and
uncertainty quantification

in geophysical modeling.
DA refers to techniques from:

statistical inference;
numerical analysis and optimization;
and dynamical systems and automatic control.

DA combines simulations from a physics-based numerical model and real-world observations of a physical process.
The output is an estimate of a posterior probability density for the numerically simulated physical state, or some statistic.
Bayesian framework: an ensemble-based forecast is a sampling procedure for the forecast-prior probability density.
Sequentially and recursively estimating the distribution for the system’s state by combining :

model forecasts, and
streaming observations

is known as filtering.

Observation-analysis-forecast cycle

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

Numerical weather prediction (NWP) uses filtering to produce weather forecasts in an:

observation
analysis
forecast

cycle.
We assume that we have an initial forecast-prior for the physics-based numerical model state.
This can represent, e.g.,

the atmosphere in terms of:

temperature,
pressure and
humidity

over Reno, six hours in the future.

The prior is generated from the simulation of the geophysical fluid dynamics.

This is described by a system of PDEs, discretized spatially over three dimensions of the atmosphere and surface topography.

Observation

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

The “true” physical process evolves in time;

eventually, the physical state reaches the time of the first forecast.
This corresponds to our real atmosphere reaching the forecast state of plus six hours.

At this time, we receive an observation of the atmosphere’s state variables.

Analysis

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

The observation comes with a likelihood function;
- we compute the likelihood of the model forecast given the observation.
Using the forecast-prior and the likelihood,

we estimate the Bayesian update of the prior to the posterior
conditioned on the observation.

The posterior density is denoted the analysis.

Forecast

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

The analysis is then used as the initialization of the next forecast.
Simulating the model forward in time, we produce a new forecast-prior.
Chaotic evolution (the butterfly effect) causes the model forecast to drift from the true state.

This is in addition to the errors in:

accurately representing the true physical process with the imperfect numerical model; and
estimating the “true” initial condition from incomplete and noisy observations.

Observation

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

When a new observation becomes available…

Analysis

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

We produce a new analysis…

Twin experiments

From: Carrassi, A. et al. Data assimilation in the geosciences: An overview of methods, issues, and perspectives. Wiley Interdisciplinary Reviews: Climate Change 9.5 (2018): e535.

And the cycle continues ad-infinitum…
Artificially generating the time-series of observations by a numerical model in parallel with the numerical model forecast, is known as a twin experiment.
In a twin experiment, a “truth-twin” is simulated to:
1. generate the artificial timeseries of observations; and
2. then to verify the accuracy of the filter’s predictions and analysis.
The “model-twin” refers to the simulated forecast used in the observation-analysis-forecast cycle.
The accuracy of the observation-analysis-forecast cycle is evaluated,

as the root-mean-square-error (RMSE) of the estimated mean state of the analysis posterior density.
This is computed using the difference of the model-twin, analysis-mean state with the known truth-twin state.

Operational constraints on filtering

An accurate representation of the true Bayesian posterior infeasible in geophysical models.
- The number of model state variables can be on the order of \( \mathbf{\mathcal{O}\left(10^{9}\right)} \);
- the number of observations on the order of \( \mathbf{\mathcal{O}\left(10^{8}\right)} \).
Therefore, the DA cycles typically estimates the first two moments of the posterior, or its mode.
- Even the storage of the covariance matrix for the model state is infeasible,
- thus empirical, low-rank covariance estimates are generated by the anomalies of the ensemble forecast.
The practical number of samples is typically on the order of \( \mathbf{\mathcal{O}\left(10^2\right)} \) due to the computational limits.
- As such, a number of approximations are introduced in filtering schemes to compensate for the severe degeneracy in the sampling.

Toy models for benchmark twin experiments

A Gaussian approximation for the:

prior,
posterior, and
likelihood

is at the basis of all operational DA schemes, due to the operational constraints.
In a linear-Gaussian model, estimating the mean or the mode is equivalent;

however, in operational DA the true distributions are often non-Gaussian and the physical dynamics are highly nonlinear.

Therefore, different filtering assumptions and approximations lead to different results in practice.
Toy models are small-scale numerical analogues of full-scale geophysical dynamics which are used in twin experiments.
Toy models:

exhibit fundamental features of large-scale dynamics;
are transparent in the model structure and assumptions;
yet are also computationally simple to simulate and sample.

Toy models are thus used for diagnostic purposes with filtering schemes;

toy models help to determine the accuracy and robustness of the various approximations and assumptions in a learning scheme under a variety of configurations.

Single layer Lorenz-96 model

Image of global atmospheric circulation.

Courtesy of: Kaidor via Wikimedia Commons (CC 3.0)

Imagine a simplified 1-dimensional model of the atmosphere around a lattitude circle.
We will suppose that this lattitude circle can be discretized into \( n \) total longitude sectors;

each sector \( i \) is represented by a single state variable \( x_i \).

The classical Lorenz-96 model is given as, \( \frac{\mathrm{d}\mathbf{x}}{\mathrm{d} t} \triangleq \mathbf{f}(\mathbf{x}), \) where,

for each state component \( i\in\{1,\cdots,n\} \), \[ \begin{align}\large{f_i(\mathbf{x}) =-x_{i-2}x_{i-1} + x_{i-1}x_{i+1} - x_i + F}.\end{align} \]

The state variables \( x_i \) have periodic boundary conditions modulo \( n \), \( x_0=x_n \), \( x_{-1}=x_{n-1} \) and \( x_{n+1}=x_{1} \).
The term \( F \) in the Lorenz-96 system is the forcing parameter that injects energy to the model.
The Lorenz-96 model approximates geophysical behavior with:

external forcing and internal dissipation with the linear terms; and
advection and conservation of energy in the quadratic terms.

Two Layer Lorenz-96 Model

From: Wilks, D. Effects of stochastic parametrizations in the Lorenz'96 system. Quarterly Journal of the Royal Meteorological Society 131.606 (2005): 389-407.

The two layer version of the Lorenz-96 system simulates coupled, ocean-atmosphere dynamics.
The inner circle represents oceanic variables evolving at a slow time scale;

these are coupled to the outter circle of atmospheric variables evolving at a fast time scale.

The two layer model is used for examining the effects of model uncertainty,

e.g., when the fast variables are not simulated accurately in the model-twin.

This is a common model to study stochastic model reduction;

instead of simulating the fast-variable dynamics, their effects on the slow variables are parameterized with a stochastic process.

This kind of approximation can be justified mathematically due to the Central Limit Theorem:

in the asymptotic separation of the time scales for the fast and slow variables,
the effect of the fast variables on the slow variables will reduce to additive, Gaussian noise which perturbs a mean-field ordinary differential equation¹.

1. Gottwald, G. et al. Stochastic climate theory. Nonlinear and Stochastic Climate Dynamics. 209–240. 2015. Cambridge University Press.

Bias in twin experiments

Two layer Lorenz-96 is used in benchmark twin experiments to study effects of model uncertainty and model reduction errors in multiscale dynamics;
- particularly, the effects of bias in the forecast prior when the physical process model generating the model-twin does not match that of the truth-twin generating the observation timeseries.
In a deterministic, biased-model setting, the numerical precision the ensemble forecast can be substantially reduced without a major deterioration of the DA cycle's (relative) predictive performance².
- The model bias overwhelms the errors that are introduced due to precision loss when the model-twin is resolved with a low order of accuracy.
However, differences in statistical properties of model forecasts of stochastic dynamical systems have been observed due to the discretization errors of certain low-order schemes.
- For example, Frank & Gottwald develop an order 2.0 Taylor scheme to correct the bias in the drift term in the Euler-Maruyama scheme in their study system³.

2. Hatfield, S. et al. Choosing the optimal numerical precision for data assimilation in the presence of model error. Journal of Advances in Modeling Earth Systems, 10, 2177–2191, 2018.
3. Frank, J. and Gottwald, G. A. A Note on Statistical Consistency of Numerical Integrators for Multiscale Dynamics. Multiscale Modeling & Simulation, 16, 1017–1033, 2018.

L96-s Model

We study how the bias of numerical integration schemes in stochastic dynamical systems also biases twin experiment statistics.

We define the L96-s model as follows, \[ \begin{align} \frac{\mathrm{d} \mathbf{x}}{\mathrm{d} t} \triangleq \mathbf{f}(\mathbf{x}) + s(t)\mathbf{I}_{n}\mathbf{W}(t), \end{align} \] where

\( \mathbf{f} \) is defined as in the single layer Lorenz-96 model
\( \mathbf{I}_n \) is the \( n\times n \) identity matrix,
\( \mathbf{W}(t) \) is an \( n \)-dimensional Wiener process; and
\( s(t):\mathbb{R}\rightarrow \mathbb{R} \) is a measurable function of (possibly) time-varying diffusion coefficients.

Both the truth-twin and model-twin are simulated with the L96-s model.

This is the analogue to a geophysical problem where the numerical model is unbiased in representing the true physics, but where the physics are intrinsically random.

This represents an idealized two layer Lorenz-96 model in which the separation of the time scales of the atmosphere and ocean is taken to infinity.

This is also a standard model configuration for filter benchmarks in the geophysical DA community.

Strong versus weak convergence

We study numerical integration schemes which converge in the strong sense;

strong convergence measures the expected path-discretization errors over all possible Wiener processes starting at an initial condition.
This is the analogue of deterministic path-discretization errors.

If \( \mathbf{x}_\mathrm{SP} \) is an exact sample path, evolving with respect to a particular Wiener process;
and \( \mathbf{x} \) is an approximation of this path, discretized at a maximum step size of \( \Delta \);
loosely, we say the approximate \( \mathbf{x} \) converges strongly to \( \mathbf{x}_\mathrm{SP} \) at order \( \gamma \) if:

there exists a \( \Delta_0 \) and a constant \( C>0 \) such that for every \( \Delta < \Delta_0 \)

then the expected path-discretization error is bounded as \[ \begin{align} \mathbb{E}\left[\left\Vert \mathbf{x}(T) - \mathbf{x}_\mathrm{SP}(T)\right\Vert\right] \leq C \Delta^\gamma \end{align} \] where these are the states, evolved from time \( 0 \) to time \( T \).

Strong versus weak convergence – continued

Weak convergence measures the error in representing some statistic of the forward distribution,
- given the evolution an initial point Dirac-delta distribution over all possible realizations of the Wiener process.
If \( g \) is a \( 2(\gamma +1) \) continuously differentiable function of at most polynomial growth;

we say that \( \mathbf{x} \) converges weakly to \( \mathbf{x}_\mathrm{SP} \) at order \( \gamma \) if for all \( \Delta< \Delta_0 \), \[ \begin{align} \left\Vert \mathbb{E}\left[ g(\mathbf{x}(T)) - g(\mathbf{x}_\mathrm{SP}(T))\right] \right\Vert \leq C \Delta^\gamma, \end{align} \] for any such statistic \( g \).

Loosely, strong convergence measures the mean of the path-discretization errors, while weak convergence can measure the error in representing the mean over all paths.

Strong versus weak convergence in twin experiments

Note: integration schemes that converge in the strong sense also converge in the weak sense,

however, weak schemes aren’t guaranteed to converge in the strong sense.

For twin experiments, there is an important distinction between these modes of convergence:

The truth-twin should be generated by a simulation which is precise in the strong sense,

here we assume we have observations of a path that is consistent with the governing dynamics ⁴.

An ensemble forecast representation of the prior needs to converge in the weak sense alone;
indeed, the forecast represents the sampling of the prior density, and we do not need to ensure the precision of each path solution if the density is accurate.

Therefore, the model-twin should be evaluated in terms of the precision in the weak sense.

4.Hansen, J. A. and Penland, C. Efficient approximate techniques for integrating stochastic differential equations. Monthly weather review. 134, 3006–3014, 2006.

Integration schemes

We study three commonly used numerical integration schemes for stochastic differential equations (SDEs):

Euler-Maruyama – a simple extension of the deterministic Euler scheme, order 1.0 strong convergence in L96-s.
4-stage Runge-Kutta – a simple extension of the determinstic 4-stage Runge-Kutta scheme⁵, order 1.0 strong convergence in L96-s
Second order Taylor – an integration scheme based on the Taylor-Stratonovich expansion⁶, order 2.0 strong convergence in L96-s.

The derivation of the order 2.0 strong Taylor scheme for the Lorenz-96 model is nontrivial and has not appeared earlier in the literature to the authors' knowledge.

Because L96-s model has:

constant or vanishing second derivatives in the deterministic component;
periodic boundary conditions; and
additive scalar noise;

we can compute this scheme efficiently, with complexity growing linearly in the system dimension \( n \).

5. Rüemelin, W. Numerical treatment of stochastic differential equations, SIAM Journal on Numerical Analysis, 19, 604–613, 1982.
6. Kloeden, P. and Platen, E. Numerical Solution of Stochastic Differential Equations, Stochastic Modelling and Applied Probability, Springer Berlin Heidelberg, page 359. 2013.

Strong convergence benchmarks

Plot of strong convergence discretization error versus step size.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. Geoscientific Model Development. In submission.

We plot the estimated strong convergence discretization error in the vertical axis versus the time-discretization step size in log-log base 10 scale.
The parameter \( s \), fixed for each panel, is the diffusion level in the L96-s model.
We note that all the orders of convergence are empirically verified here.

However, the constant \( C \) in the bound \( C\Delta^\gamma \) has a large impact on the overall strong discretization error.

Particularly, we see that in the low diffusion regime, the order 1.0 Runge-Kutta scheme has an overall discretization error comparable to the order 2.0 Taylor scheme.
The performance of the Runge-Kutta scheme varies significantly, however, between low and high diffusion.

While the order of convergence remains the same in all diffusion regimes, the difference in performance is reflected in the constant \( C \).

Weak convergence

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. Geoscientific Model Development. In submission.

The effect of the constant is even more pronounced in the weak convergence metric.
The order 1.0 Runge-Kutta scheme outperforms the order 2.0 Taylor scheme in low diffusion.
The overall weak discretization error varies on an order of magnitude between low and high diffusion levels.
High precision in the strong sense across all diffusion regimes makes the Taylor scheme a natural choice for generating the truth-twin.
The Runge-Kutta scheme is appealing to generate the ensemble forecast of the model-twin, due to its accuracy in the weak sense.

However, because the weak convergence of the Runge-Kutta scheme varies by an order of magnitude between diffusion levels, we use the Taylor scheme as a benchmark in the following experiments.
We study how the ensemble forecast statistics of the Euler-Maruyama and the Runge-Kutta schemes differ from those generated by the Taylor scheme.

Summary of bias in ensemble forecast

Plot ensemble forecast bias versus time.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. Geoscientific Model Development. In submission.

For lack of time, we summarize our benchmark of ensemble forecast statistics here.
We confirm the robustness of Runge-Kutta in generating ensemble forecast statistics in the L96-s model.
For a fine discretization of \( \Delta=10^{-3} \) there is virtually no difference in the ensemble forecast statistics of the Runge-Kutta scheme from the benchmark.
Even for a coarse time step of \( \Delta=10^{-2} \), has unbiased spread compared with Taylor using a step size of \( \Delta=10^{-3} \).
While divergence of the Runge-Kutta trajectories from Taylor trajectories occurs earlier with a coarse time step,

divergence of their ensemble means matches asymptotic divergence of the finely resolved Runge-Kutta using time-step of \( \Delta=10^{-3} \).

Summary of bias in ensemble forecast

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. Geoscientific Model Development. In submission.

Euler-Maruyama introduces systematic biases into ensemble forecast statistics.
The asymptotic ensemble spread is systematically, artifically larger.
There is extreme divergence of trajectories causing the ensemble mean to depart from the benchmark Taylor scheme;
the divergence is over an order of magnitude difference compared with the Runge-Kutta scheme.
This is concerning as Euler-Maruyama is commonly used to simulate systems of SDEs for twin experiments.
This indicates the biasing of ensemble forecast statistics,

however, this doesn’t yet demonstrate the effect of this bias on a filter twin experiment.

This is demonstrated in the following experiements.

Taylor benchmark configuration – filter statistics

Plot of benchmark Taylor configuration filtering statistics.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. Geoscientific Model Development. In submission.

We plot the twin-experiment RMSE and ensemble spread for an ensemble Kalman filter (EnKF);
both the truth-twin and model-twin are generated with the order 2.0 Taylor scheme.
Both twins use step size \( \Delta=10^{-3} \).
Vertical axis: diffusion level \( s \).
Horizontal axis: variance of the error in the observations given to the ensemble Kalman filter.

We vary the method of ensemble integration between Runge-Kutta and Euler-Maruyama, with a step size between \( \Delta\in\left\{10^{-3},10^{-2}\right\} \).

Filter benchmark – Runge-Kutta model-twin, coarse discretization, Taylor truth-twin, fine discretization

Plot of coarse Runge-Kutta filter simulation versus Taylor benchmark.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. Geoscientific Model Development. In submission.

The difference between generating the ensemble with the Runge-Kutta with a step size \( \Delta= 10^{-3} \) and benchmark Taylor is negligible.

Filtering statistics for these two configurations differ on the order of \( 10^{-6} \).

Increasing the step size of Runge-Kutta to \( \Delta=10^{-2} \), we introduce small errors,

however, the residuals are unstructured in sign or magnitude.

The mean of the differences in the RMSE is approximately \( 8\times 10^{-5} \) though with a standard deviation of approximately \( 10^{-3} \) (on the order of the weak-discretization error).
This indicates that the observed differences in the filtering statistics can be attributed to random numerical errors, which are themselves unbiased.

Filter benchmark – Euler-Maruyama model-twin, fine discretization, Taylor truth-twin, fine discretization

Plot of of coarse Runge-Kutta filter simulation versus Taylor benchmark.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. Geoscientific Model Development. In submission.

Generating the ensemble forecast with Euler-Maruyama with time step \( \Delta=10^{-3} \), we find structure in the residuals.
In low diffusion, the RMSE using Euler-Maruyama suffers from biases observed in the ensemble forecast.
The spread suffers from the artificial inflation in the forecast.
However, in high diffusion, the RMSE and spread of the ensemble are nearly identical using either Taylor or Euler-Maruyama.

The level at which the bias impacts the twin experiment statistics depends strongly on the observation uncertainty, and most especially on the model uncertainty.
We demonstrate the effect of increasing the step size of Euler-Maruyama to \( \Delta=10^{-2} \) in generating the ensemble.

Filter benchmark – Euler-Maruyama model-twin, coarse discretization, Taylor truth-twin, fine discretization

Plot of coarse Euler-Maruyama filter simulation versus Taylor benchmark.

From: Grudzien et al. On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiments. Geoscientific Model Development. In submission.

The structure of the residuals remains largely the same for the coarse discretization, but the magnitude changes dramatically.
In low diffusion the forecast discretization error forces filter divergence.
In high diffusion, \( s\geq 0.5 \), the difference of the RMSE with the benchmark Taylor is on the order of \( 10^{-2} \).

We consistently see the bias in the ensemble spread, present once again.

Conclusions

Typically, Lorenz-96 is simulated with a 4-stage Runge-Kutta scheme with step size up to \( 0.5 \) in deterministic settings.
- We show that the same standards cannot apply to twin experiments in the commonly used form of the L96-s.
We must distinguish between strong and weak convergence, and its impact on the truth-twin and the model-twin respectively.
Euler-Maruyama is a commonly used integration scheme for SDEs, but we find that it introduces strong, systematic bias into twin experiments when the step size is greater than or equal to \( \Delta=10^{-2} \).
- Effects of the bias depends strongly on the observation and especially model uncertainty in a twin experiment;
- however, the bias is sufficient to cause filter divergence in weak diffusion.
We find that the 4-stage Runge-Kutta scheme is a statistically robust solver, without the systematic biases encountered in the Euler-Maruyama scheme.
- For small step sizes \( \Delta \) on \( \mathcal{O}\left(10^{-3}\right) \), we find virtually no difference between the Runge-Kutta scheme and the higher order Taylor scheme.
- A step size of \( \Delta=10^{-2} \) introduces additional discretization error, but the error doesn't significantly influence the RMSE or spread of the EnKF.

Conclusions – continued

Due to the computational constraints in running twin experiments, a step size of \( \Delta=10^{-3} \) is often unreasonable for generating the truth-twin or model-twin.
A step size on \( \mathcal{O}\left(10^{-2}\right) \) is what can be afforded in a typical benchmark study.
This excludes the use of Euler-Maruyama from statistically robust twin experiments entirely.
We demonstrate in our work that an overall discretization error can be bounded by \( 10^{-3} \), using:

4-stage Runge-Kutta with step size \( \Delta=10^{-2} \) to generate the model-twin; and
order 2.0 Taylor with step size \( \Delta=5\times 10^{-3} \) to generate the truth-twin.

This forms a practical compromise, which our diagnostics demonstrate does not bias the outcomes of the filtering statistics.
We provide a computationally efficient framework for statistically robust twin experiments in the L96-s model.
We provide a novel derivation of ther order 2.0 Taylor scheme for the L96-s model.

This derivation has not previously appeared in the literature to the authors' knowledge, and we moreover provide benchmarks on the convergence of this and other schemes.

As an open question:

Can using weak integration schemes (which may not converge to any sample path) reduce the computational burden of ensemble-based forecasts in geophysical models?

If the goal of the forecast is to converge in distribution alone, this may be a viable alternative to traditional geophysical modeling paradigms.